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-j  Finite  element  algorithms  are  presented  £or  the 
approximate  solution  of  the  streamfunction-vortlcity 
equations  of  steady  incompressible  viscous  flows.  Both 
the  linear  Stokes  and  the  nonlinear  Navier-Stokes 
equations  are  considered.  The  methods  discussed  re¬ 
quire  low  continuity  finite  element  spaces  and  do  not 
require  any  artificial  specification  of  the  vorcicity 
at  solid  boundaries.  Particular  attention  is  paid  to 
methods  for  multiply  connected  domains  and  to 
theoretical  and  computational  estimates  for  the 
accuracy  of  the  algorithms.  Brief  consideration  is 
also  given  to  three  dimensional  problems,  to  exterior 
problems,  and  to  the  recovery  of  (he  pressure  field. 
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where  p  is  the  constant  density,  we  have  that,  for 

i  ■  1, . . . ,m, 

0  •  i-  Vp-T_  do  *  (vAu-^)*Vii  +  £)  «t  da. 


Since  fm-T  »  3<j/3n  whenever  div  u  ■  0  and 
ii-Vu-T_  ■  uii *n  +  S[(u.u)/2]/3T,  we  then  have  that 

|  (v  o)£*n  +  £-t)dc  *  0  for  1  ■  1 . m. 


Thus  (1,2, 5, 6,  and  8)  are  the  governing  equations  and 

The  stationary  Navier-Stokes  equations,  written  in  side  conditions  which  are  to  determine  the  functions  V p 

of  the  streamfunction  i {/  and  the  vortlcity  u>,  and  a  and  the  constants  a^,  1  *  l,...,m. 


terms  of  the  streamfunction  iji  and  the  vortlcity  0), 
are  given  by 

Ail)  -  -0)  in  (2  (1 


Weak  Formulation 


We  define  the  function  spaces 


-vAu  +  <§■  S  -  S  * curl  i  in  n  (2) 

where  in  (1)  and  (2)  v  is  the  kinematic  viscosity,  £ 
the  body  force,  and  ft  is  a  bounded  region  in  R  .  If 
n  is  multiply  connected,  we  denote  by  T0  its 
exterior  boundary  and  by  Tj,  i  “  l,...,m,  the  re¬ 
maining  parts  of  the  boundary.  See  Figure  1. 

Suppose  that  at  these  boundaries  the  velocity  u 
isspecified,  i.e.,  u-£  on  r-u?.nri-  order 

for  a  streamfunction  to  exist,  we  must  have  that  [1] 


Hr(fl)  -  {*  e  L2(f2); 


e  L2(fi) ,  a,B,p,r  e  Z+, 


da  -  0  for  i 


P.  P  “  0, 


i.e.,  there  is  no  nec  mass  flow  through  any  of  the 
boundary  pieces.  It  is  well  known  that  the  stream- 
function  is  unique  up  to  an  additive  constant,  and  we 
fix  its  value  by  specifying  it  to  be  zero  at  an 
arbitrary  point  xQ  on  ro.  Now,  let  q  denote  a 
function  such  that 

1^-  *  g-n  on  T  and  q(x  )  *  0.  (4) 

Then  the  boundary  conditions  for  the  system  (1)  and  (2) 
are  given  by 

ij>  »  q  on  i(>  »  on  1^,  i  -  l,...,m  (5) 


U  “  -4*T  on  T,  (6) 

where  a^,  i  *  l,...,m,  are  constants  to  be  determined 
as  part  of  the  solution.  These  constants  are  fixed  by 
the  requirement  that  che  pressure  be  single  valued, 
i.e.,  the  change  in  the  pressure  p  around  any  closed 
curve  surrounding  any  of  i  •  l,...,m,  is  zero. 

Since  the  momentum  equation  stipulates  that 

—  7p  »  -u-7u  +  \i£u  +  f  (7) 

0  y  -  -  “  “ 
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H1  (12)  -  {0  E  H1  (12) ;  <J>  —  0  on  D, 

O 

H^((2)  -  {*  e  H1  (12) ;  »  ■  0  on  y  *  ■  c±  on  1^, 
i  •  l,...,m,  c^  arbitrary), 

V1,4^)  -  t<t>  e  L4(J2);  ||.|£  e  L4 (J2) } , 

V  ■H1ffl)n«1,4(ll),  V-H1(I1)i1«1,4(I1), 

a  a 

and  the  set 

Sa  -  {<»>  e  H1  ((2)  n  W2’4  (12) ;  $  -  q  on  To; 

$  -  q+c  on  f  ,  i  »  1 . .  c  arbitrary). 

i  i  1  1 

r*  2 

In  addition  we  will  make  use  of  the  spaces  H  (D 
consisting  of  traces  of  functions  belonging  to  Hr(I2) 
for  r  E  Z f,  and  also  the  dual  spaces  with  negative 
indices  of  the  above  defined  spaces.  For  details  con¬ 
cerning  these  spaces,  see  [1,2]. 

The  weak  formulation  of  the  streamfunction- 
vortlcity  problem  which  we  will  utilize  is  as  follows. 

Given  f  £  [L2((2)]2  and  4  £  (H  2(0]2  satisfying 
(3),  we  seek  ij)  £  Sa  and  u  e  V  such  that 

f  WC  d!2  -  |  curl  ili-curl  I  df2  •  Lg/X  ds  (9) 

i  h~ 

for  all  LEV, 
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The  Linear  Case 


-v  [  curl  wcurl  0  dfi  +  [  |^)dn  (10) 

j  -  j  -\v  3x  :X  3v 

a  a 

*  -  f  *curl  $  dil  for  all  0  e  V 

1 - 

where  curl  -  (3/3y,-3/3x) .  This  weak  formulation  is 
the  most  practical  special  case  of  a  more  general 
formulation  found  In  [1,31.  Also,  this  formulation  was 
used  in  [4]  to  successfully  compute  high  Reynolds 
number  driven  cavity  flows  on  nonuniform,  relatively 
coarse,  grids.  We  note  at  the  outset  that  the  only 
boundary  conditions  explicitly  imposed  on  Che  functions 
appearing  in  (9)  and  (10)  are  those  corresponding  to 
(5).  In  particular,  (6)  is  a  natural  boundary 
condition  and  no  boundary  conditions  on  u  are  imposed 
on  boundaries  where  (5)  and  (6)  apply.  Furthermore, 
the  constraints  (8)  are  also  natural  to  the  formulation 
(9-10) .  To  elucidate  these  points  and  to  show  the 
connection  between  (9-10)  and  (1,2, 5, 6,  and  8),  let  us 
proceed  formally  and  perform  appropriate  integrations 
by  parts  in  (9-10).  We  are  then  led  to 


We  will  also  be  interested  in  the  linear  Stokes 
flow  case,  and  it  will  be  advantageous  to  sometimes 
consider  this  case  separately.  For  this  case,  one 
simply  omits  all  terms  arising  from  the  nonlinear  con¬ 
vection  terms  in  (2).  Thus  (2)  is  replaced  by 
Cthi  -  -curl  £  where  we  have  absorbed  the  constant  v 
into  {_.  Also,  the  weak  formulation  (9-10)  is  replaced 
by  the  following.  We  seek  e  H^(I2)  and 
0  E  {0  e  Hl(0);  0  ■  q  on  ro;  0  ■  q+Cj  on  Tj,  i  *  l,...,m, 
cA  arbitrary}  such  that  (9)  and 

f  curl  tii -curl  $  d(i  «  [  f  -curl  $  dfl  (13) 

n  n 

for  all  0  e  H^(fi) 

are  satisfied.  Note  chat  it  is  no  longer  required  for 
0,0  E  W^,4((2)  since  this  inclusion  results  from  the 
need  to  make  the  nonlinear  term  in  (10)  well  defined. 
Also,  (8)  is  replaced  by 


and 


j  (u»A0K  dft 
a 


A 

l3n 


+  4*£)C  do, 


(ID 


+  i*i)d0 


0  for  i  ■  1,. 


.  ,m. 


(14) 


|  (vAw~0  ~  +  |^  ||  +  curl  f)0  dfi  (12) 

SI 

“  |  (V  -  (i)0  +  f.*X)$  do. 

r 

Since  we  may  choose  C  to  alternately  vanish  and  not 
vanish  on  the  boundary  I",  we  recover,  from  (11),  (1) 
and  (6).  Also  (2)  Is  recovered  from  (12)  by  choosing 
0  to  vanish  on  the  boundary  T.  Since  0  Is  required 
to  vanish  on  ro,  the  integral  on  the  right  hand  side 
of  (12)  vanTsbes  on  ro.  By  alternately  choosing  0 
to  vanish  on  all  1  -  l,...,m,  except  on  one,  say 

Ti,  on  which  0  Is  required  to  be  a  constant,  that 
same  Integral  yields,  in  view  of  (4)  and  (5),  that  (8) 
Is  satisfied.  Then  Indeed  the  side  conditions  (6)  and 
(8)  are  natural  to  the  formulation  (9-10),  and  the  only 
essential  boundary  conditions  are  those  on  0  Itself, 
i.e.,  (5).  Again  we  note  that  no  boundary  condition 
on  oi  need  be  Imposed. 

In  addition  to  the  ease  with  which  the  boundary 
conditions  are  satisfied,  the  formulation  (9-10)  also 
allows  for  the  use  of  low  continuity,  i.e.,  merely 
continuous  over  3,  function  spaces  and  thus  greatly 
simplifies  its  finite  element  discretization.  Also, 
as  a  practical  observation,  we  note  that  the  function 
q  is  required  only  on  the  boundaries  T.,  i  *  0,..,,m, 
and  thus  usually  may  be  easily  computed  from  (4) 


IX  -  Discretization 


The  discretization  of  (9-10),  or  of  (9,13)  in  the 
linear  case,  proceeds  in  the  usual  manner,  with  the 
only  difficulty  resulting  from  the  Inhomogeneity  q  in 
the  essential  boundary  condition  (5) .  If  the  boundary 
F  consists  of  polygons,  then  one  merely  chooses  a 
finite  dimensional  subspace  Vh  C  V,  and  then  define 
the  space 

Vh  •  {0  e  V1*;  0  -  0  on  r  ;  0  *  cj  on  F, , 
a  oil 

1  "  l,...,m,  c^1  arbitrary) 

and  the  set 

sj|  “  {0  e  V*1 ;  0  •  q*1  on  0  ■  q^+cj  on  T^, 
i  »  l,...,m,  c^  arbitrary) 


where  q  E  V 
since  one  need 
may  choose  qh 
of  q  in  Vh 
0h  £  SjJ  and  Co 


is  an  approximation  to  q.  For  example 
define  q“  only  along  boundaries,  we 
on  T  to  be  the  boundary  interpolant 
restricted  to  T.  Then  one  seeks 
h  E  V*1  such  that 


separately  on  each  of  these  parts  of  the  boundary. 

So  far  we  have  only  considered  boundary  conditions 
whict^  correspond  to  a  specification  of  the  velocity  at 
all  boundaries  and  we  will  continue  to  focus  on  this 
case  In  the  sequel.  However,  by  examining  (11-12),  we 
see  that  many  other  kinds  of  boundary  conditions  can 
be  also  implemented.  For  example,  suppose  f  »0  and 
that  is  «  3u/3n  ■  0  on  part  of  the  boundary,  say  T0. 
Such  a  boundary  condition  is  useful  in  matching  a 
viscous  flow  to  an  external  lnvlscid  flow.  In  this 
case  we  retain  (9-10)  with  f.  *  0  but  now 
a,;  E  H^O)  are  required  to  vanish  on  T0  and  0,  0 
are  not  specified  there.  With  u  *  0  on  f0  and  0 
arbitrary  there,  we  see  that  the  right  hand  side  of 
(12)  yields  3u/3n  •  0  as  a  natural  boundary 
condition. 


f  dfl  -  |  curl  0h  •  curl  C*1  dfl 


*  do  for  all  c  V*1, 


,  h  .  ^  f  ,  h'-Ah  30^  3  ih  30h, 

-v  curl  0)  *curl  $  (-r - r- - -7 - r-— )dTi 

-  -  j  3X  cX  0y 


-  f -curl  0h  d(. 


(15) 


(16) 


Curved  boundaries  may  be  treated  by  isoparametric 
finite  elements  and  similar  techniques.  The  linear 
Stokes  equations  can  be  discretized  in  the  same  manner 
with  Vh  C  Hi  Cl) . 

The  discrete  weak  formulation  (15-16)  is  equiva¬ 
lent  to  a  system  of  nonlinear  algebraic  equations.  The 
latter  is  arrived  at  by  choosing  a  basis  {$j). 

j  *  1.....J,  for  V*1.  Suppose  we  order  these  basis 
functions  so  chat 


*  0  on  T  for  J  -  1,...,N 


*  -  0  on  r-Tj  for  j  -  N+M^+l . »4«1  and 

for  1  -  0, . . . ,m, 

where  •  0  and  N+M^  “  J.  Thus  the  first  N 

basis  functions  correspond  to  Interior  degrees  of 
freedom  or  to  degrees  of  freedom  on  r  which  are  not 
associated  with  function  values  and  the  basis  functions 
indexed  by  ]  •  N+Mj_j+1, . . .  correspond  to  the 

degrees  of  freedom  associated  with  function  values  on 
I*.  For  example,  if  Vh  is  a  Lagrangian  finite 
element  space,  Chen  for  j  -  1,...,N  correspond 

to  interior  nodes  while  for  j  ■  . . .  .N+M^ 

correspond  to  nodes  on  T^.  Thus,  in  this  example,  H 
is  the  number  of  nodes  in  the  interior  of  (1  and 
are  the  number  of  nodes  on  Ti  for 
i  «  0,...,m.  He  note  that  in  a  practical  imple¬ 
mentation  one  would  not  choose  to  number  the  basis 
functions  according  to  (17).  He  merely  choose  that 
numbering  scheme  in  order  to  simplify  the  exposition 
which  follows. 

Having  chosen  the  numbering  scheme  (17),  we  may 
then  substitute  in  (15-16) 


<*)h  -  l  (*) 
1-1  J  3 


.  N  J 

V  -  l  +  [  <l(£, >M*) 

j-1  3  3  J-N+l  ^  J 


+  l  l  v 

1-1  1  l"1’4*!-]*1  J 


where  xj  denotes  Che  coordinates  of  the  j-th  node  and 
where  we  have  approximated  q  on  T  by  its  boundary 
interpolant.  He  may  also  choose,  in  (15-16), 

-  ik(x)  for  k  -  1,...,J  and  -  $^(x)  for 
k  -  1,...,S  and  (x)  where  the  sum  ranges 

over  )  •  ft+M^j+l , . . . , for  i  -  l,...,m.  Thus 
(15-16)  represent  (fW-J+a)  nonlinear  algebraic  equations 
for  the  (N+J+m)  unknowns  <u, ,  J  -  1,...,J,  ♦j  , 

J  -  1,...,N,  and  a^,  1  -  l,...,m.  Belov  we  remark 

on  how  to  solve  this  discrete  system. 

For  Che  case  of  the  linear  Stokes  problem,  the 
weak  formulation  (15-16),  where  now  the  nonlinear  terms 
in  (16)  are  not  present,  leads,  in  an  analogous  manner 
to  a  linear  system  of  (N+J+m)  algebraic  equations. 

Nonlinear  Solver 

The  nonlinear  system  of  equations  resulting  from 
(15-16)  nay  be  solved  in  a  variety  of  ways.  He  are 
articuiarly  Interested  in  preserving  the  feature  that 
ic  boundary  conditions  on  the  vorticity  be  imposed  at 
boundaries  for  which  the  velocity  is  known.  This 


precludes  the  use  of  algorithms  which  Iterate  between 
(15)  and  (16).  Even  so,  one  nay  choose  to  imbed  (15-16) 
into  a  real  or  pseudo-time  dependent  problem  and  then 
discretize  the  tine  derivatives,  or  one  may  choose  to 
use  a  nonlinear  equation  solver  such  as  Newton's  method, 
a  quasi-Newton  method,  or  even  a  fixed  point  method. 
Because  of  its  robustness  and  ease  of  programming,  we 
use  Newton's  method,  although  qussi-Newton  methods, 
when  applicable,  may  be  more  efficient.  He  here  give 
the  Newton  algorithm  which  is  most  easily  described  in 
terms  of  the  weak  formulation.  To  further  simplify  the 
presentation  we  introduce  the  bilinear  forms 


A(w,C)  -  f  u;  dn  for  u,s  e  V 
1 
n 


B('l'.C)  -  -f  curl  '1-curl  ;  dfi  for  £  V,  (20) 

J 

£1 


the  trlllnear  form 


C(<a.lM)  -  |  u(f£f£  -  |||*>dS)  for  E  V>  (21) 


and  the  linear  functionals 


F($)  -  -  f-curl  $  d£2  and  G($)  -  f  ifcg’T  d£)  for  $  c  V. 

I  l 


Then,  (15-16)  take  the  form 


A(o)h,Ch)  +  B(*h,Ch)  -  GUh)  for  all  Ch  E  Vh,  (23) 


vB($h,oh)  +  C(0)h,i|jh,(j>h)  -  F(<t>h)  for  all  <t>h  E  vj\  (24) 


Newton's  method  for  (23-24)  is  then  given  as  follows. 
Given  u°  E  Vh  and  ip°  £  S^,  we  generate  the  sequence 
{wn,<J/D}  for  n  >  1  by  solving  the  linear  system 


A(un,Zh)  +  B«/\;h)  -  G(Ch)  for  all  ;h  E  Vh 


vB($h,iun)  +  C(ajn.tfrn-1,«J>h)  +  C(wn-1,Hin,*h)  (26) 

-  F($h)  +  C(un~1,vn_1,$h)  for  all  <J>h  E 

Due  to  the  relatively  small  attraction  ball  of  Newton’s 
method,  it  is  convenient  to  start  with  the  following 
simple  iteration  method  which,  at  least  at  low  Reynolds 
numbers,  is  globally  convergent.  Given  tc°  E  Vh  and 
E  va’  w*  g®11*18**  the  sequence  for  n  1 

by  solving  the  linear  system 


A(un.i:h)  +  BWn,;h) 


G(;h)  for  all  ;h  c  Vh  (27) 


vB(ifh,u)n)  +  C(UB-1.*n,:h)  -  F(0h)  for  all  ih  E  Vh.  (28) 


Thus  the  composite  algorithm  consists  of  doing  a  few 
steps,  usually  one  or  two,  of  (27-28)  in  order  to  get 
into  the  neighborhood  of  the  solution  of  (23-24),  and 


then  switching  to  Newton's  method  (25-26)  in  order  to 
sore  quickly  hose  in  on  that  solution.  Noce  chat  the 
simple  iteration  algorithm  (27-28),  and  therefore  the 
composite  algorithm  as  well,  do  not  need  the  initial 
guess  to  satisfy  the  boundary  conditions.  In 

particular,  we  may  choose  %j°  «  <u°  •  0. 

By  methods  similar  to  chose  used  in  [1,5,6]  for  the 
primitive  variable,  i.e.,  velocity-pressure,  formulation 
it  can  be  shown  that  the  Newton  iterates  defined  by  (25- 
26)  converge  quadracically  to  a  solution  of  (23-24)  for 
a  sufficiently  close  initial  guess.  Further  remarks 
concerning  the  solution  of  the  discrete  nonlinear 
system  (23-24)  are  made  below  when  computational  results 
are  discussed.  We  noce  that  either  of  the  methods 
(25-26)  or  (27-28)  require  the  solution  of  a  sequence 
of  linear  algebraic  systems. 

Two  Mechods  for  Multiply  Connected  Domains 

As  was  previously  discussed,  the  main  difficulty 
which  multiply  connected  domains  present  is  that  the 
streamfunction  may  be  arbitrarily  specified  on  only  one 
part  of  the  boundary.  Thus,  above,  we  have  set  i]/  and 


at  one  point  on 


but  these  are  determined 


only  up  to  the  unknown  constants  a^  and  aj, 
respectively,  on  the  remaining  parts  Tit  i  ■  l,...,m, 
of  the  boundary.  In  the  discretization  method  repre¬ 
sented  by  (15-16),  the  constants  a^  are  determined  as 
part  of  the  solution  process.  However,  that  method 
requires  the  use  of 


l  ♦.(*),  1  -  l....,m  (29) 

j-N+Mi_l+l  J 

as  basis  functions  associated  vith  the  a!j\  i  ■  l,...,m. 
In  general,  these  basis  functions  are  only  "semi-local" 
in  Che  sense  that  they  couple  all  the  points  on  a 
boundary  T..  In  particular,  in  the  discrete  set  of 
equations,  a!j  is  coupled  with  all  the  unknowns 
associated  with  elements  whose  closures  intersect  vith 
Tj.  The  result  of  this  is  that  the  bandwidth  of  the 
linear  systems  encountered  will  be  greatly  increased, 
resulting  in  both  a  computer  storage  and  time  penalty 
over  simply  connected  domain  problems.  Some  of  this 
penalty  may  be  mitigated  by  employing  a  numbering 
schema  resulting  in  banded- bordered  matrices  wherein 
Che  onerous  couplings  are  found  near  the  bottom  and  the 
right  of  the  matrices  encountered.  This,  of  course, 
requires  the  development  of  special  programs  to  solve 
the  linear  systems. 

An  alternative  Co  using  basis  functions  such  as 
(29)  is  co  guess  the  value  of  the  constants  a^, 

i  -  1 . m,  appearing  in  (5).  Then  one  may  solve  for 

and  ;  from  (1,2,5  and  6).  However,  in  general, 

(8)  will  not  be  satisfied.  At  this  point  one  may 
change  the  guesses  for  the  a^'s,  repeating  the  process 
until  convergence  is  achieved,  i.e.,  until  (8)  is 
satisfied.  In  more  detail,  we  proceed  as  follows. 

Given  guesses  a^,  i  •  1 . a,  k  e  Z +,  we  define 

c  sk  *n<*  c  vh  where 


^•(M  H1  (0)  nn1'*®;  9  *  q  on 


q+a*  on  rit  i  -  1 . m}, 


to  be  the  solution  of  (15-16)  where  nov  (16)  holds  for 
all  t  Vg  •  Vh  ~  H^(f.).  In  this  case  is  given 

by  (IS)  with  aj  replaced  by  the  known  numbers  aj 
and  thus  the  discrete  system  (15-16)  contains  (N+J) 
equations  for  the  (N+J)  unknowns  a,,  j  *  1,...,J,  and 


'lj,  j  *  Furthermore,  these  equations  and  un¬ 

knowns  may  be  ordered  in  such  a  wav  that  the  bandwidth 
of  the  discrete  system  is  no  larger  than  that  obtain¬ 
able  in  an  analogous  discretization  of  a  problem  posed 
in  the  simply  connected  domain  bounded  by  ro.  Having 
obtained  from  a£,  i  •  l,...,m,  we  may  substitute 

into  (8).  In  practice,  we  would  evaluate  the  Integrals 
appearing  in  (8)  by  a  numerical  quadrature  formula. 

Let  us  denote  such  an  approximation  to  each  equation  in 
(8)  by  I^fa*1) ,  i  -  l,...,m,  which  in  general  will  not 
be  small.  Of  course,  each  Ij(*)  is  a  nonlinear 
function  of  a^  »  (aj,...,aJj)  through  the  discrete 
vorticity  tu|j.  We  now  need  an  algorithm  to  generate 
new  guesses  £*c+^,  from  which  and  <ojj+1  may  be 

computed  and  for  which  I^(a^+^)  is  closer  to  zero 
than  was  I^(£^).  The  most  practical  way  to  accomplish 
this  is  to  use  a  secant  type  method.  For  example,  if 
m  ■  1,  i.e..  Cl  is  doubly  connected,  we  may  then 
define 


1  1  _  ,  k. 

T  ,  k.  ,  k-1.  W' 
1l(al)*1l(al  > 


For  m  >  1,  some  generalized  secant  method  in  R  such 
as  the  Wolfe  secant  method  can  be  used.  See  [7]  for 
details  about  such  methods.  Note  that  the  use  of 
formulas  such  as  (30)  requires  two  starting  guesses  a0 
and  £*.  We  also  note  that  in  general  it  is  not 
practical  to  use  Newton's  method  to  update  £k  since 
evaluating  the  Jacobian  of  1(a)  -  (1^,...,]^,)  requires 
the  solution  of  linear  problems  of  the  same  size  as 
(15-16)  for  each  pair  3u)/3 aj,  at. 

Thus,  at  the  price  of  solving  a  sequence  of 
simpler  problems,  one  may  avoid  the  bandwidth  problems 
engendered  by  the  use  of  basis  functions  such  as  (29). 
Whether  or  not  the  alternate  method  is  more  efficient 
depends,  of  course,  on  how  many  of  the  simpler 
problems  must  be  solved  as  well  as  on  how  much  cheaper 
it  is  to  solve  the  simpler  problems.  Boch  of  these 
factors,  for  a  particular  geometry,  will  depend  on  the 
number  of  degrees  of  freedom,  i.e.,  the  grldsize,  and 
the  Reynolds  number.  Therefore  it  is  clear  that  the 
performance  of  iterations  such  as  (30)  play  a  central 
role  in  the  overall  efficiency  of  the  alternate  method 
of  treating  multiply  connected  domains.  Unfortunately, 
the  sequence  £k  generated  by  secant  type  methods  are 
not  guaranteed~to  converge  for  arbitrary  values  of  the 
initial  guesses,  especially  at  high  values  of  the 
Reynolds  number.  We  will  return  to  this  point  when  we 
discuss  some  computational  results  below. 

Multiply  Connected  Domains  in  the  Linear  Case 

The  situation  for  the  alternate  method  discussed 
above  is  much  simpler  in  the  case  of  the  linear  Stokes 
equations.  In  this  case,  I(£)  is  a  linear  function 
of  the  a^'s  and  thus,  given  any  £°  and  £*,  a  secant 
method  such  as  (30)  will  yield  that  £^  is  the  exact 
desired  value  £^.  Furthermore,  and  are 

linear  combinations  of  uj)  and  tjj.  k  ■  0,1,  and  thus 
the  former  may  be  obtained  from  the  latter  without 
solving  another  linear  problem  of  the  type  (15-16) 
where  in  (16)  the  nonlinear  terms  are  omitted. 

These  observations  may  be  used  to  define  an  even 
simpler  algorithm  for  the  linear  case  which  is  equiva¬ 
lent  to  the  above  alternate  method.  First  we  solve 
(m+1)  discrete  problems  corresponding  to  the  continuous 
problems 


and  -vl-^  *  i  in  0, 


\  -  qk  on  r#.  'I*k  -  qk+ck  on  ri,  i  -  1, 


vk  k 

-r-  -  -I  on  r 


where  k  «  0,..,,m  and  f°  -  f_,  q°  •  q,  g°  -  4, 
c°  “0,  £k  -  q*  m  4k  .  0,  and  cj  ■  5^  for  k  >_  1 

and  i  •  1 . m.  All  of  these  problems  involve  the 

same  left  hand  side,  i.e.,  the  same  coefficient  matrix 
in  the  discrete  problem,  and  thus  may  be  solved 
simultaneously  by  setting  up  (nri-1)  right  hand  sides. 
Now  consider  the  combination 


$  -  <P  +  l  a^ij),  and  5  -  «  +  J  a.  i|/k  (32) 
k-1  k”l 


where  a^,  k  _  l,...,m,  are  constants  to  be  determined, 
then,  since  the  problems  (31)  are  linear,  we  have  that 


u-uj  .  We  will  consider  the  case  of  q  -  0  and  g  •  0 
and  of  simply  connected  domains.  Inhomogeneous 
problems  can  be  treated  by  techniques  similar  to  those 
used  in  [8]  for  the  primitive  variable  formulation  of 
Navier-Stokes  equations  and  insofar  as  the  accuracy  of 
the  approximate  solution,  problems  posed  on  multiply 
connected  domains  should  behave  in  a  manner  similar  to 
those  posed  on  simply  connected  domains.  Furthermore, 
we  will  focus  mostly  on  low  order,  low  continuity 
finite  element  spaces,  especially  continuous  piecerwlse 
linear  polynomial  spaces. 

Most  of  the  results  available  are  concerned  with 
the  linear  Stokes  problems  in  a  simply  connected  domain 
with  homogeneous  boundary  data.  Indeed,  finite  element 
discretizations  of  the  particular  weak  form  (9,13)  are 
considered  in  [9-14].  All  except  for  [13]  and  [14] 
consider  only  the  case  of  piecewise  polynomial  spaces 
of  degree  two  or  higher.  The  analyses  of  [14]  Improves 
on  the  previous  works,  and  it  is  the  results  of  that 
work  which  we  summarize  here. 

To  begin  vith,  we  define  the  space  of  weakly 
ha  '.sonic  functions 


{;  E  H1  (12)  ;  VC*V<J>  dft 


for  all  ♦  E  H1  (!))}. 

0 


and  -v&u  •  curl  f  la  0 


i>  *  q  on  r  , 


q+a1  on  r±,  i  -  1, 


We  also  define  the  norms 


j  an  |  u 

l!*ILl  "  a"d  "*"*  ’  *eH 


£--4-1  -  r. 

Thus,  for  any  values  of  the  constants  at,  $  and  u 
satisfy  the  given  Stokes  problem  (1,5,6  and  13).  We 
fix  the  aj's  by  requiring  that  (14)  be  satisfied. 
Indeed,  denoting  each  of  the  integrals  in  (14)  by 
I1(4i),  we  see  that  since  these  are  linear  in  u 


0  -  I1(m)  -  IjO^)  +  I  *icI1(“k^  for  1 


Thus  (33)  is  a  linear  system  of  m  equations  for  the 
n  constants  a^,  k  ■  l,...,m.  Having  found  these 
constants,  then  (32)  yields  the  solution  of  the  given 
Stokes  problem. 

Of  course,  the  use  of  the  continuous  problems  (31) 
was  symbolic;  in  reality  one  solves  the  corresponding 
discrete  problems  and  approximates  the  integrals 
w  by  numerical  quadrature.  In  summary,  we  see 
that  in  the  linear  case,  we  may  solve  multiply 
connected  domain  problems  by  solving  a  single  "simple" 
matrix  problem,  with  multiple  right  hand  sides, 
evaluating  the  (m+l)m  integrals  l^^),  k  *  0,...,m, 
1  *  1, ...,a,  solving  the  a  *  m  linear  system  (33), 
an ^dliheu  taking  the  linear  combinations  (32).  On  the 
otmh  hand,  one  may  solve  the  same  problem  by  using 
basis  functions  such  as  (29).  In  this  case  one  need 
solve  a  single  linear  system  with  a  single  right  hand 
side  to  obtain  the  solution.  However,  as  noted  above, 
this  system  vlll  surely  be  more  complicated  than  that 
for  the  method  (31-33),  and  it  seems  that,  at  least 
fo*jpedarate  values  of  m,  the  latter  technique  is 
preferable. 


Pjpoderati 
» (arable. 


where  ||  •  |[  for  r  E  Z  denotes  a  norm  on  H  (ft). 
Note  that  tne  first  of  these  is  not  the  norm  on  H~^(ft), 
the  dual  space  of  hJ(0),  but  we  adopt  this  notation 
for  simplicity.  A  basic  result  is  that  the  two 
norms  of  (34)  are  equivalent  on  H.  (The  proof  of  this 
and  the  other  statements  concerning  the  linear  case  may 
be  found  in  [14].) 

Now,  recall  chat  our  approximations  are  based  on 
choosing  a  subspace  V*1  C  H^-(ft)  which  we  will  assume 
to  be  a  continuous  piecewise  polynomial  space,  and  then 
letting  V*J  »  Vh  O  Hg(n).  We  may  then  define  the  space 
of  discrete  weakly  harmonic  functions 


curl  C^'curl  dfi 


0  for  all  d>h  e  A 


In  general  W*1  C  H.  In  fact,  the  lack  of  this  in¬ 
clusion  results  in  the  non-optimality  of  the  approxi¬ 
mations  to  u  and  also  causes  the  main  difficulty  in 
the  analyses.  Having  defined  Hh,  we  may  also  define 
Che  norm 

[  Ah  dll 


In  general  the  norma  I!  *  II  * , h  and  II  *  ere  not 
equivalent  on  H*1.  However,  for  any  E  H"  and 

E  >  0  we  have  that  for  some  constants  Ci 


We  now  turn  to  a  brief  discussion  of  the  available 
theoretical  estimates  of  the  differences  and 


-.\v.vV.* 


We  next  define  the  differences 


h  ,h  ~h  *?h  h  h  _  .  „ 

e^  -  v  -v  .  e1  -  y-1  ,  e2  ■  u  -R^  and  e^  -  oi-ILu  (37) 


.T.h  .  „h 


where  ;Ji  £  Vq  is  arbitrary  and 


u  c  yb 


satisfies 


|  curl  (lo-R^a)  •  curl  dfl  «  0  for  all  $h  e  v|\ 

3  (38) 

Also,  we  noce  that  froa  (9,13)  with  £  •  0  and  q  •  0 
and  from  its  discrete  analogue, 

[  curl  ei?-curl  dfl  •  0  for  all  $h  e  V*1  (39) 

J  -  2  -  o 

a 

j  (e^;h  +  curl  e^ »curl  Ch)dfl  •  f  (e.,lh  curl  e^  *curl  S^dR 

n  a 

for  all  (‘tv11,  (40) 


He  then  obtain  that  for  any  e  >  0  and  for  some 
constants  Cj,  and  Cj, 


||u-(oh|lo  <_  ^l!u-Rh<JH0  +  ^ 


NcuglQMOlL  <  2  inf  Hcurl(»-»  )ll 
*h  J 
*  eVo 


+  C3ll«**-Rhw||o  +  C4  +  C5  h1/2_eE^ 


llourl  eJ;'o  <  Jle^lLj  +  !i  !!_1  +  I'curl  e2  !?o  (46) 


-  Cll|e2!'*,h  +  C2  **1/2_e ‘l*2^o  +  lle2!!-l 


+  curl  e. 


Then  (42)  follows  from  (37),  (45),  and  the  triangle 
inequality,  and  (41)  follows  from  (37),  (44-46),  and 
the  triangle  inequality. 

If  Hh  C  H,  then  Eg  -  e£  -  0.  Thus  (41)  and  (42) 
yield  that 

J|a>-tuh !| o  <_  2||w-Rhw||o 


||curl(p»-’i>h)  II  <  2  inf  ||curl(i|i-$h)  II  +  C  ||«-IL'j|] 

0  “  ;h  „h  o  3  h  o 

<ii  ev 

o 
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which  are  optimal  estimates  for  the  L  -norm  of  the 
vorticity  and  of  the  velocity  field.  Unfortunately, 
the  inclusion  H*1  C  H  implies  that  v"CH  (D),  i.e., 
Vh  consists  of  continuously  differentiable  piecewise 
polynomials.  In  this  case  we  might  as  well  have  dis¬ 
cretized,  using  V*1  C  H2(R),  the  fourth  order  stream- 
function  formulation  of  the  Stokes  equations. 

Of  more  interest  to  us  here  is  the  case  of 
V*1  d  H2(R),  e.g.,  Vh  consisting  of  merely  continuous 
piecewise  polynomials.  Examining  (41)  and  (42)  yields 
that  except  for  the  terms  involving  Eg  and  eJ,  the 
terms  on  the  right  hand  sides  are  either  approximation 
theoretic,  or,  due  to  (38),  can  be  approximated  in 
terms  of  approximation  theoretic  results  [14],  Thus  it 
only  remains  to  estimate,  for  specific  choices  of  Vh, 
the  terms  Eg  and  Eg.  For  example,  if  Vh  consists 
of  piecewise  linear  polynomials  and  if 
¥  e  H2(R)  n  Hq(0),  then  for  any  e  >  0  and  j  *  0,1, 

It  can  be  shown  that 


C6+J  h 


J  +I-e 


curKii-ji  )  curl  5  dfl 


Er  «  inf  sup 

MM 


for  j  •  0,1. 


Then  (41)  and  (42)  yield,  for  linear  finite  element 
spaces. 


|o>-whI|o  -  0(h1/2)  and  ||curlOH|ih)  ||q  -  0(h1_e) .  (47) 


To  prove  (41)  and  (42),  first  note  that  by  (39), 
e^  £  H*1  and  from  (40)  with  Ch  e  fr*  one  easily 
obtains  that 


i*2n*.h  i  ■M-i +  •up„ 


>2”o  i  Mo  + 

'hctr 


curl  e^»curl  5  dR 


curl  e^ -curl 


We  note  that  (47)  shows  that  the  velocity  field 
u  *  curl  is  optimally  approximated  and  that  there  is 
a  3/2’ s  loss  in  the  exponent  of  h  in  the  estimate 
for  4>. 

For  the  nonlinear  case,  there  are  much  fewer  error 
estimates  available  in  the  literature.  In  fact,  for 
the  particular  method  considered  here,  such  estimates 
can  only  be  found  in  [1].  Furthermore,  these  hold  only 
for  polynomial  spaces  of  degree  two  or  higher,  and 
yield  suboptimal  results.  However,  by  combining  the 
nonlinear  analysis  of  [1]  with  the  improved  analysis, 
for  the  linear  problem,  found  in  [14],  results  such  as 
those  above  can  be  reasonably  expected  to  hold  for  the 
nonlinear  problem  as  well. 

Ill  -  Remarks 


Also  (36)  and  (40)  yield  that 


Infinite  Domain  Problems  -  So  far  our  considerations 
have  focused  on  problems  which  may  be  posed  on  bounded 


WWW 
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domains.  However,  the  streamfunction-vorticltv  formu¬ 
lation  may  be  especially  useful  for  problems  posed  on 
exterior  domains  since  for  many  such  problems  the 
vorticity  decays  exponentially  with  the  distance  from 
the  origin  while  the  velocity  and  pressure  only  decay 
algebraically.  Thus  applying  far  field  conditions  on 
the  vorticity  can  lead  to  substantially  smaller  compu¬ 
tational  regions  than  that  needed  for  a  computation  of 
the  same  accuracy  using  such  conditions  on  Che  velocity. 
We  have  also  indicated  in  Section  I  how  easily  imple¬ 
mented  boundary  conditions  on  the  vorticity  may  be 
useful  in  matching  to  external  inviscid  flows  in 
boundary  layers  and  wake  type  calculations. 

Finite  Difference  Methods  -  It  would  be  remiss  not  to 
point  out  that  finite  difference  methods  may  also  be 
implemented  in  such  a  manner  that  no  artificial  boundary 
condition  on  the  vorticity  is  imposed  at  solid 
boundaries.  Indeed,  the  key  to  avoiding  such  a  specifi¬ 
cation  is  not  the  use  of  finite  element  methods. 

Rather  it  is  the  willingness  to  solve  Che  discrete 
system  (15-16)  as  one  coupled  set  of  equations  as 
opposed  to  iterating  between  (15)  and  (16).  Thus,  in 
a  finite  difference  method,  one  could  discretize  (1) 
and  (6)  while  making  use  of  (5)  as  well,  yielding  more 
equations  than  unknowns  determining  ^  as  a  function 
of  u  and  the  boundary  data.  One  also  discretizes 
(2)  without  imposing  any  boundary  condition  on  m, 
yielding  less  equations  than  unknowns  to  "determine" 
u.  In  this  manner  one  may  not  solve  the  discrete 
versions  of  (1,S,6)  for  <|>  since  that  system  is  over- 
constrained,  and  also  one  may  not  uniquely  solve  the 
discrece  version  of  (2)  for  <u  since  that  system  is 
underconstrained.  However,  the  above  type  of  discreti¬ 
zations  of  (1,2,5  and  6),  taken  together,  may  be  solved 
simultaneously  for  i|»  and  m.  We  note  that  similar 
observations  about  the  relation  of  the  discrete  systems 
resulting  from  (1,5,6)  and  (2)  hold  in  the  finite 
alement  case. 

The  method  of  treating  multiply  connected  domains 
wherein  one  uses  iterations  such  as  (30)  can  clearly  be 
used  in  connection  with  finite  difference  methods  where 
again  one  must  approximately  evaluate  the  constraint 
(8)  or  (1A)  in  order  to  update  the  guessed  values  of 
ip  at  the  boundaries.  Also  in  a  straightforward 
manner,  one  may  implement  the  other  method.  Specifi¬ 
cally,  one  can  leave  the  constants  at  in  (5)  as  un¬ 
knowns,  and  then  add  a  discretization  of  (8)  or  (1A)  in 
order  to  close  the  system  of  equations.  However,  we 
note  that  in  the  finite  element  case  (8)  or  (1A),  as 
well  as  (6),  are  natural  boundary  conditions,  and  thus 
are  more  easily  satisfied.  The  same  remarks  concerning 
the  relative  efficiency,  e.g.,  bandwidth  size,  of  the 
two  methods  which  held  for  the  finite  element  case  also 
hold  in  the  finite  difference  case. 


that,  by  integrating  by  parts,  that  (A8)  nay  be  viewed 
formally  as  a  discretization  of  the  problem 


Ap  -  o  div(-u/Vu  +  vAu  +  f) 


^  “  p(-u-Vu  +  vAu  +  f)  -n 


Both  of  these  may  be  obtained  directly  from  (7)  by 
respectively  taking  the  divergence  of  (7)  and  the  inner 
product  of  (7)  with  £.  Unfortunately,  the  right  hand 
side  of  (A8)  is  not  defined  when  one  uses  merely 
continuous  finite  element  spaces  for  iph  and  M*1.  The 
problem  is  not  with  the  viscous  terms  since,  from  the 
definition  of  u  ,  we  have  that 


f  Auh-Vph  d  ft 


f  „  3 
P 


(A>P  )  do 


where  the  last  equality  is  only  approximate.  The 
Integral  on  the  right  is  well  defined,  even  for  low 
continuity  spaces.  The  problematical  term  is  the 
convection  term  in  (A8)  which  is  not  defined  for  merely 
continuous  spaces.  However,  for  the  linear  case  or 
the  nonlinear  case  with  finite  element  spaces  which  are 
at  least  continuously  differentiable,  one  may  use  (A8), 
with  the  replacement  (A9),  to  solve  for  the  pressure. 

We  now  present  a  method  for  recovering  the 
pressure  which  works  in  all  cases.  First,  (7)  may  be 
expressed  in  the  form 


VH  -  p[(_^)u>  +  £+vAu]  where  H  -  (p+i  pu*u),  (50) 

l.e.,  H  is  the  total  pressure  head.  We  then  consider 
the  problem  of  seeking  H  E  H-1  (0)  D  W^’"(ft)  such  that 


j|  VH-VH  dft  "  p  (-Mu -curl  H  +  f^VH  +  v7H-Au)dft 

“  p  j  (-Mu-curl  H  +  f/VH)dft  -  pv  j  H  da 


for  all  H  e  (ft)  O  W4,1(ft) 


Recovery  of  the  Primitive  Variables  -  The  computation 
of  the  velocity  field  from  the  streamf unction  is  a 
simple  matter  since  we  may  define  u11  •  (3iJ/h/3y,-3iiih/3x) . 
The  recovery  of  the  pressure  field  is  not  so  straight¬ 
forward,  especially  for  low  continuity,  e.g.,  merely 
continuous,  finite  element  spaces.  Formally,  one 
would  like  to  use  (7)  in  the  following  way.  From  17), 
we  define  a  discrete  pressure  ph  E  V“  C  HMft)  n  Vr ’"(ft) 
such  that 

f  7ph-Vfch  dT.  -  p  [  (-uh-Vuh+vVuh+f).V£h  dft  (A8) 

h  ft 

for  all  £h  £  V11 


which  by  additional  integration  by  parts,  is  formally 
equivalent  to 


AH  »  p  curl (uu)  +  p  div  £  in  ft 


If  ■  p  (.1*2. _  v  If  +  u2.*i)  °°  r  • 


These,  using  div  u  *  0,  may  be  derived  directly  from 
(50).  Thus  we  are  led  to  define  an  approximation 
Hh  E  the  same  finite  element  space  used  for  the 
vorticity  approximation,  by  the  solution  of  the 
following  problem.  We  seek  H*1  £  V*1  such  Chat 


vhare  uh  in  the  right  hand  side  is  obtained  from  the 
ahead v~eeapu ted  approximate  streamfunctlon  .  Note 


f  VHh-?Hh  d.l  -  3  f  <£-uih7iih)»7Hh  dr. 


-  PM  Hh  da  for  all  Hh  e  Vh 


where  now  the  right  hand  side  Is  well  defined  for  Che 
already  computed  values  of  e  V*1  and  t|<"  £  Sc¬ 
once  the  approximate  total  head  H*1  is  computed  from 
(52),  the  pressure  Is  easily  recovered  from 
ph  .  yh  .  ((^h)2  +  (ily) 2 ] /2 .  We  noce  that  the  solution 
H“  of  (52)  corresponds  to  an  approximate  solution  of 
the  Netasann  problem  (51).  Thus  H*1  and  therefore  p*1 
are  determined  only  up  to  an  additive  constant,  which 
is  to  be  expected. 

Three  Dimensional  Problems  -  The  streamf unction- 
vortlcity  formulation  of  the  Navier-Stokes  equations 
has  not  achieved  the  same  interest  or  success  in  three 
dimensional  settings  as  it  has  in  two  dimensions. 
However,  recently  there  has  been  increasing  attention 
devoted  to  such  problems.  See,  e.g.,  [15,16]  for 
finite  difference  approximations  of  three  dimensional 
problems. 

Since  div  ii  ■  0,  we  have  that  necessarily 
u  -  curl  for  some  vector  valued  function 
variously  called  che  'Vector  streaafunction"  or  the 
"vector  velocity  potential".  Of  course,  the  vorticlty 
is  defined  to  be  u  «  curl  u.  Then,  the  Navier-Stokes 
equations,  in  terms  of  the  streamf unction  £  and  ’ 
vorticlty  oj,  are  given  by 


curl  curl  i  “ 

v  curl  curl  £  “  curl  -  u*V(curl  jfO 

«  curl (w  x  curl  £). 


At  a  boundary  where  ij  is  specified,  one  would  now 
specify  curl.<£. 

As  a  consequence  of  their  definition,  it  follows 
Chat  <[>  can  be  determined  only  up  to  the  gradient  of 
an  arbitrary  scalar  function,  and  that  div  u  -  0. 

There  are  various  ways  In  which  these  facts  have  been 
used  to  simplify  the  formulation.  For  example,  the 
arbitrariness  In  £  can  be  pinned  dawn  by  requiring 
that  div  <j>  •  0.  This  method,  popular  In  electro¬ 
magnetic  problems  where  it  is  called  the  "Coulomb 
gauge”,  has  been  used  in  [15]  in  connection  with  a 
finite  difference  solution  of  vortex  flows  in  all  of 
R^,  i.e.,  a  problem  with  no  boundaries.  In  problems 
with  solid  boundaries,  it  becomes  difficult  to  enforce 
boundary  conditions.  The  obvious  disadvantage  of  this 
method  is  that  there  are  six  unknown  scalar  fields. 

The  main  advantage  of  this  method  is  that  (53)  and 
div  »  0  imply  Chat  -A£  “  w,  so  that  together  with 
(54),  the  governing  equations  may  be  viewed  as  coupled 
#oisson  equations  for  the  six  scalar  fields  constituting 
Che  components  of  £  and  u.  Although  more  fields  ar  *. 
req^red  than  in  the  primitive  variable  formulation 
and^e  many  as  in  the  veloclty-vortlcity  formulation, 
the  resulting  streamfunction  vorticlty  equations  are 
presumably  easier  Co  solve.  Furthermore,  the  finite 
elamanp  techniques  discussed  in  this  paper  for  the 
plane  flow  setting  extend  in  a  straightforward  manner 
Co  this  three  dimensional  method. 

Another  way  of  fixing  the  streamfunction  is  to 
set  one  component,  say  the  component  'P,  in  the 
x-dlrectlon,  to  zero.  The  obvious  advantage  of  this 
method,  which  was  used  successfully  in  [16]  to  compute 
compressible  flows,  la  that  we  have  only  two  unknown 
streamfunction  fields.  In  addition,  as  Is  Indicated  In 
[16J[,  this  method  can  handle  solid  boundaries  better 


than  the  first  technique  and  nay  be  especially  useful 
for  flows  which  are  perturbations  of  two  dimensional  or 
axially  symmetric  flows.  In  [16]  a  third  technique  is 
also  discussed,  namely  letting  it  «  7*  *  Vc  which  also 
involves  only  two  scalar  fields.  However,  this  choice 
introduces  additional  nonlinearities  into  the  problem 
which,  in  practice,  is  not  desirable. 

IV  -  Computational  Examples 

Accuracy  -  We  first  examine,  through  computational 
examples,  the  accuracy  of  the  particular  finite  element 
methods  discussed  in  Section  II.  The  context  of  these 
examples  is  smooth  solutions  of  the  linear  Stokes 
equations  posed  on  the  unit  square,  i.e.,  a  simply 
connected  domain.  However,  we  note  that  the  algorithms 
have  been  successfully  used  to  compute  solutions  of  the 
nonlinear  Navier-Stokes  equations.  For  example,  in  [4] 
these  methods  are  used,  in  conjunction  with  a  reduced 
basis/continuation  technique,  to  compute  accurate 
approximations  of  driven  cavity  flows  at  Reynolds 
numbers  up  to  10,000  using  relatively  coarse  nonuniform 
meshes. 

Some  of  the  computational  results  for  the  Stokes 
equations  are  summarized  in  Tables  1  and  2  which 
respectively  deal  with  piecewise  linear  and  piecewise 
quadratic  finite  element  spaces  for  both  i|>  and  <u. 
Each  table  gives  the  exact  solutions  for  lji  and  (ii. 

In  some  of  the  cases  of  Table  2  this  required  the 
introduction  of  an  lnhomogenelty  in  (1),  i.e.,  we  have 
that  Aijj  +  to  “  a  for  some  function  a.  Other  infor¬ 
mation  contained  in  the  tables  is  whether  or  not  the 
boundary  conditions  on  \ p  and  3i|)/3n  are  homogeneous, 
i.e.,  whether  or  not  q  and  g  in  (5)  and  (6)  vanish, 
and  also  whether  a  uniform  grid  or  a  graded  nonuniform 
grid  was  used  in  the  calculations.  Finally,  each  table 
contains  the  computed  rates  of  convergence  of  the 
finite  element  approximations  to  and  to  as 
measured  in  the  L^(C1)  and  H^(fl)  norms,  i.e., 
respectively  the  mean  square  errors  in  the  function 
values  and  in  the  derivatives.  These  rates  were 
computed  by  comparing  errors  on  different  grids. 

The  most  obvious  trend  in  these  results  is  that 
the  streamfunction  \J>  and  its  derivatives  are  always 
optimally  approximated.  On  the  other  hand,  there  is 
in  general  a  loss  of  accuracy  in  the  vorticlty  approxi¬ 
mation.  In  particular,  there  seems  to  be  a  loss  of  one 
power  of  h,  the  grid  size,  in  the  piecewise  linear 
case,  and  a  loss  of  the  3/2  power  of  h  in  the 
quadratic  case.  Thus  the  theoretical  results  of 
Section  II  seem  to  be  sharp  in  both  cases  for  the 
streamfunction  and  also  for  the  vorticlty  in  the 
piecewise  quadratic  case.  We  also  note  that  in  the 
piecewise  linear  case  we  always  obtain  optimal  approxi¬ 
mations  to  p  and  0)  whenever  all  boundary  conditions 
are  homogeneous,  e.g.,  see  1-4  in  Table  1.  This  was 
not  the  case  for  the  piecewise  quadratic  case,  e.g., 
see  7  in  Table  2.  Finally,  we  have  also  found  that 
whenever  the  approximation  to  p  is  exact,  i.e.,  'p 
belongs  to  the  approximating  space,  then  the  approxi¬ 
mation  to  ai  is  optimal.  See,  e.g.,  13  and  14  in 
Table  2.  This  result  can  be  gleaned  from  (41)  since 
in  this  case  •  0. 

Multiply  Connected  Domains  -  We  also  report  on  some 
preliminary  computations  for  problems  posed  on  multiply 
connected  domains,  namely  with  the  view  of  comparing 
the  two  methods  of  treating  such  problems  and  of 
examining  their  behavior  as  parameters,  e.g.,  the  grid 
size  or  Reynolds  number,  are  varied.  The  domain  .1 
is  the  unit  square,  i.e.,  T0  is  the  boundary  of  that 
square,  from  which  we  have  removed  a  rectangle,  i.e., 

Tj  is  the  boundary  of  a  rectangle  contained  within 
the  unit  square.  Thus,  we  deal  with  a  doublv  connected 
domain.  In  all  che  computations  we  will  have  that  q 
on  T|  vanishes  so  that  (5)  requires  that  ;  ■  a  * 
unknown  constant  on  chat  part  of  the  boundary. 


These  preliminary  computational  results  are 
summarized  in  Table  3.  All  computations  were  performed 
using  uniform  grids  of  size  h.  Results  for  three 
values  of  the  Reynolds  numbers  Re  *  1/v  are  given  and 
for  both  methods  of  treating  multiply  connected  domains, 

l.e.,  using  iterations  such  as  (30)  to  update  guessed 
values  of  the  streamfunctlon  on  Che  boundaries  or  using 
"semi-local"  basis  functions  such  as  (29)  to  directly 
compute  the  approximate  solution.  The  three  problems 
considered  can  be  characterized  by  the  boundary  value 
of  i)  at  the  left  and  right  boundarr.es  and  the 
position  of  the  hole.  Specifically,  we  consider  the 
following  problems: 


X.  iHO.y)  -4Kl,y)  -2y3-3y2 


3.  i^(0, y)  -  <1/(1, y)  •  y 


(1/4, 1/2)  x  (1/4, 3/4) 


fij  -  (1/4.1/2)  x  (1/4, 1/2) 
S\  -  (2/5, 3/5)  x  (2/5, 3/5) 


where  2^  is  the  region  bounded  by  The  remaining 

boundary  conditions  on  ^  are  constant  values  at 
y  •  0  and  1  and  3<p/dn  *  0  on  all  boundaries.  We 
note  that  the  solution  of  Problem  1,  due  to  symmetries, 
should  have  'ji  ■  a  «  -1/2  on  F^.  We  have  also 
computed  with  the  boundary  condition  of  Problem  1  or  2 
with  the  concentric  hole  %  •  (1/4, 3/4)  x  (1/4, 3/4), 
and  due  to  the  high  symmetry  of  this  configuration,  the 
exact  value  of  a  ■  -1/2  was  computed  for  all  values 
of  Re  and  h  and  for  both  methods.  We  also  note 
Chat  for  the  Re  -  0  cases,  l.e.,  for  the  linear 
Stokes  equations,  the  results  using  the  iterative 
technique  always  converged  in  the  expected  one  step. 
Also,  for  the  iterative  method,  the  initial  guesses  for 
the  secant  iteration  (30)  were  a°  -  1  and  a1  -  2/3 
in  all  cases. 

From  these  preliminary  results  It  seems  that  at  a 
given  value  of  h  che  direct  method  (B)  yields  better 
accuracy  than  the  iterative  method  (A) .  Also  compu¬ 
tations  at  higher  Reynolds  numbers  sometimes  resulted 
in  the  lack  of  convergence  of  the  iterative  method  for 
the  above  initial  values.  This  is  indicative  of  the 
possible  convergence  problems  of  the  iterative  method. 
Further  and  more  detailed  computational  results  will 
be  reported  on  elsewhere. 
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Figure  1.  A  multiply  connected  domain. 
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Table  1.  Computational  resales  using  piecewise  linear  polynomials.  H  *  homogeneous, 
I  “  inhomogeneous ,  U  *  uniform,  N  *  nonuniform.  Exact  -j  ■  -lv. 


Exact  solution 


Boundary  conditions 

i   3iii/3n 


Rates  of  convergence 
in  L2  :a  in  L2  li  in  H2 


,2  2 
.  sin  ttx  sin  Try 

H 

H 

U 

2 
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.  same 

H 

H 

N 

2 

2 

.  x2y2(x-l)2(y-l)2 

H 

U 

N 

2 

2 

.  1.+  3. 

H 

H 

U 

2 

? 

.  cos  Ty 

I 

H 

U 

2 

i 

.  x2(x-l)2 

I 

H 

u 

2 

i 

Table  2.  Computational  results  using  piecewise  quadratic  polynomials.  H  ■  homogeneous, 
I  -  inhomogeneous.  All  examples  use  a  uniform  grid. 


Table  3.  Computational  results  for  doubly  connected  region  problems  using  piecewise 
linear  functions;  A  ■  using  iterative  updating  of  ^  on  boundaries, 

B  »  using  "semi-local"  basis  functions. 
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